#!/usr/bin/env python
#
#########################################################################
#
#Module: v.to.r3strata based on r.plane
#PURPOSE: Creates r3 volume map from a vector point map
#This module is adapted from the r.plane module that shipped with GRASS7.0
#July 19, 2013
#Several differences
##First the input is pulled from a vector points layer rather than user query
##This generates a series 
##Second the output cells are r3 volumes rather than r2, in order to get a cell value as well as xyz coordinates
#
#
#
#by Tim Bailey, Humboldt State Univesity
#
#
#########################################################################


import sys 
import os 
import math 
import string 
import grass.script as grass


def main():

    name = options ['output']#need to add naming procedure for generating multiple r2 raster horizons 
    type = options ['type']
    dip = float(options['dip'])
    az = float(options['azimuth'])
    ea = float(options['easting'])
    no = float(options['northing'])
    el = float(options['elevation'])
    cv = int(options['cellvalue'])      #int for category cell values, r.plane produces, cv is assigned from the point and the entire r2 has the same value
    doi = int(options['distance_of_influence']      # not yet implemented, intended to limit inference using a distance
             
    dsm = 
            ###raster surface elevation
    ed = float(options['euclidean distance']

    reg = grass.region()

    ###test input values
    if abs(dip)>= 90:
                      grass.fatal(_("dip must be between -90 and 90 degrees."))
    if az < 0 or az > 360:
                      grass.fatal(_("azimuth must be between 0 and 360"))
    

    
    ###r.plane algorithm
    az_r = math.radians(-az)        ####azimuth is counterclockwise whereas standard notation for geology at least is clockwise for positive azimuth
    sinaz = math.sin(az_r)
    cosaz = math.cos(az_r)

    dip_r = math.radians(-dip)
    tandip = math.tan(dip_r)

    kx = sinaz * tandip
    ky = cosaz * tandip
    kz = el - ea * sinaz * tandip - no *cosaz * tandip
#   ed = ((kx - ea)^2 + (ky - no)^2 + (kz - el)^2)^.5)
               
    ()

######
####enforcing nulls cells greater than surface elevation.  If kz > dsm(z) then null
# 

    grass.mapcalc("$name = $type($round(x() * $kx + y() * $ky + $kz +cv))"
                  name = name, type = dtype, round = round, kx = kx, ky = ky, kz = kz  cv = cv)
#not sure if this is going to successfully genereate
               
#   grass.mapcalc("$name = $type(               if $kz > dsm, null else
                  
#####enforcing limits of distance to doi.  If doi > ed then null

    
    grass.run_command('r.support', map = name, history = '')
    grass.raster_history(name)

    grass.message(_("Done."))
    t = string.Template("Raster map <$name> generated by v.to.r3strata

    



    
    









    

    
    

    
